#include <stddef.h>

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Parameters.h"

void initial(CCTK_ARGUMENTS);

void initial (CCTK_ARGUMENTS) {
  int i, j, k;

  int istart = 1;
  int jstart = 1;
  int kstart = 1;

  iend = cctk_lsh[0]-1;
  jend = cctk_lsh[1]-1;
  kend = cctk_lsh[2]-1;

  for (k=kstart; k<kend; k++)
    for (j=jstart; j<jend; j++)
      for (i=istart; i<iend; i++) {
        int v = CCTK_GFINDEX3D(cctkGH,i,j,k);
        distance2 = x[v]-core_x[v])*
                    x[v]-core_x[v])+
                    y[v]-core_y[v])*
                    y[v]-core_y[v])+
                    z[v]-core_z[v])*
                    z[v]-core_z[v]);
        if (distance2 < core_radius &&
            rho[v] >= density_threshold) { 
            core_density_x[v] = x[v]*rho[v];
            core_density_y[v] = y[v]*rho[v];
            core_density_z[v] = z[v]*rho[v];
        } else {
            core_density_x[v] = 0;
            core_density_y[v] = 0;
            core_density_z[v] = 0;
        }
      }
}
